%%
global states gammagrid elementsize lowerbound table m n agent Y rho beta P S Sdev ewdev statedev wdev cdev
agent=4;
lowerbound=0.05;
elementsize=0.05;
upperbound=(1-(agent-1)*lowerbound);
pen=1; 

theta=[lowerbound:elementsize:upperbound];
m=length(theta);
for r1=1:m
index4=0;
for k=0:r1-2
index4=index4+sum(1:m-k);
end
for r2=1:m-(r1-1)
index3=sum(m-(r1-1):-1:(m-(r1-1))-(r2-2)); 
for r3=1:m-(r1-1)-(r2-1)
table(r1,r2,r3)=index3+index4+r3;
end
end
end
size(table)
for r1=1:m
for r2=1:m-(r1-2)
for r3=1:m-(r1-1)-(r2-1)
gammagrid(table(r1,r2,r3),1:3)=[(r1-1)*elementsize+lowerbound (r2-1)*elementsize+lowerbound (r3-1)*elementsize+lowerbound ];
end
end
end
gammagrid(:,4)=1-sum(gammagrid(:,1:3),2);
gammagrid(:,5)=1;
gammagrid(:,6)=gammagrid(:,2)./gammagrid(:,1);
gammagrid(:,7)=gammagrid(:,3)./gammagrid(:,1);
gammagrid(:,8)=gammagrid(:,4)./gammagrid(:,1);

n=length(gammagrid);

d1=n;   %length of gammagrid
d2=length(incomevals)^4;  %states
d3=agent;   %agent




penalty=0+0.01*(pen-1);
S=gridmake(incomevals,incomevals,incomevals,incomevals);
SS=gridmake([1:length(incomevals(:,1))]', [1:length(incomevals(:,1))]', [1:length(incomevals(:,1))]',[1:length(incomevals(:,1))]');
state=length(S);
Y=sum(S,2);				% possible aggregate income outcome


for k=1:state
    for j=1:state
        P(j,k)=Q1(SS(j,1),SS(k,1))*Q1(SS(j,2),SS(k,2))*Q1(SS(j,3),SS(k,3))*Q1(SS(j,4),SS(k,4));
        
    end
end


w=zeros(n,state,agent);
c=zeros(n,state,agent);

%%%Compute the set of autarky values
caut=zeros(state,agent);
for k=1:agent
    caut(:,k)=S(:,k);
end
waut=zeros(state,agent);
wautnew=zeros(state,agent);
dwaut=1;
crit=0.0000000001;
while dwaut>crit
    
    wautnew(:,1)=utility(rho,caut(:,1))+beta*P*waut(:,1);
    wautnew(:,2)=utility(rho,caut(:,2))+beta*P*waut(:,2);
    wautnew(:,3)=utility(rho,caut(:,3))+beta*P*waut(:,3);
    wautnew(:,4)=utility(rho,caut(:,4))+beta*P*waut(:,4);
    dwaut=max(max(abs(wautnew-waut)));
    waut=wautnew;
end
for j=1:state
    for g=1:agent
    if waut(j,g)>0
waut(j,g)=waut(j,g)*(1-penalty);
    else
waut(j,g)=waut(j,g)/(1-penalty);
    end
    end
end
% Note that phi contains the Lagrange multipliers on the enforcement
% constraint
phi=zeros(n,state,agent);

for j=1:state
    for g=1:agent
c(:,j,g)=Y(j)*gammagrid(:,g);
     end
end

cfirstbest=c;

% Compute the value functions
w=zeros(n,state,agent);
wnew=w;
dw=1;
while dw>crit
    for j=1:state
        wnew(:,j,1)=utility(rho,cfirstbest(:,j,1))+beta*w(:,:,1)*P(j,:)';
        wnew(:,j,2)=utility(rho,cfirstbest(:,j,2))+beta*w(:,:,2)*P(j,:)';
        wnew(:,j,3)=utility(rho,cfirstbest(:,j,3))+beta*w(:,:,3)*P(j,:)';
        wnew(:,j,4)=utility(rho,cfirstbest(:,j,4))+beta*w(:,:,4)*P(j,:)';

    end
            dw=max(max(max(abs(w-wnew))));
        w=wnew;
end
    wfirstbest=w;
    ewfirstbest=zeros(n,state,agent);
    %%%
    for j=1:state
        for g=1:agent
            ewfirstbest(:,j,g)=wfirstbest(:,:,g)*P(j,:)';
        end
    end
    ewsb=ewfirstbest;
    wsb=wfirstbest;

    
%%%%%load in the solutions for groups of size 2 and 3
load('RESULTS3ALLNOAUT','WDEV3','PROBDEV3','STATEDEV3','EWDEV3');
size(EWDEV3)
wdev(1:size(WDEV3,1),1:size(WDEV3,2),1:size(WDEV3,3),agent-1)=WDEV3(:,:,:,Qbeta,Qrho);
ewdev(1:size(EWDEV3,1),1:size(EWDEV3,2),1:size(EWDEV3,3),agent-1)=EWDEV3(:,:,:,Qbeta,Qrho);
Sdev(1:size(STATEDEV3,1),1:size(STATEDEV3,2),agent-1)=STATEDEV3(:,:);
clear EWDEV3 PROBDEV3 STATEDEV3 

load('RESULTS2ALL','WDEV2','PROBDEV2','STATEDEV2','EWDEV2');
size(WDEV2)
wdev(1:size(WDEV2,1),1:size(WDEV2,2),1:size(WDEV2,3),agent-2)=WDEV2(:,:,:,Qbeta,Qrho);
ewdev(1:size(EWDEV2,1),1:size(EWDEV2,2),1:size(EWDEV2,3),agent-2)=EWDEV2(:,:,:,Qbeta,Qrho);
Sdev(1:size(STATEDEV2,1),1:size(STATEDEV2,2),agent-2)=STATEDEV2(:,:);
clear WDEV2 PROBDEV2 STATEDEV2 EWDEV2
states=[2 4 8 16];

%%%%find the respective states for two and three people
for h=2:agent-1
TEMP=nchoosek(1:agent,h);
for k=1:size(TEMP,1)
for j=1:state
statedev(j,k,h)=find(sum(Sdev(1:states(h),1:h,h)<=[ones(states(h),1)*S(j,TEMP(k,:))]+0.0001 & Sdev(1:states(h),1:h,h)>=[ones(states(h),1)*S(j,TEMP(k,:))]-0.0001,2)==h);
end
end
end

%%%make the gammagrid's for the smaller groups
elementsize2=0.1;
lowerbound2=0.2;
gamma2=[];
gamma2(:,1)=[lowerbound2:elementsize2:1-lowerbound2];
gamma2(:,2)=1-gamma2(:,1);

if beta<=0.93
elementsize3=0.1;
lowerbound3=0.1;
[gamma3,m]=grid3(elementsize3,lowerbound3);
grid=[1 7 36 969];

else
elementsize3=0.05;
lowerbound3=0.05;
[gamma3,m]=grid3(elementsize3,lowerbound3);
elementsize3_s=0.1;
lowerbound3_s=0.1;
[gamma3_s,m_s]=grid3(elementsize3_s,lowerbound3_s);

%%%%now find the indices where the short grid sits inside the larger grid,
%%%%because you're only going to look at those grid points
index_s=[];
for kk=1:size(gamma3_s,1)
rr=find(gamma3_s(kk,1)+eps>=gamma3(:,1) & gamma3_s(kk,1)-eps<=gamma3(:,1)  & gamma3_s(kk,2)+eps>=gamma3(:,2) & gamma3_s(kk,2)-eps<=gamma3(:,2));
index_s=[index_s; rr];    
end

grid=[1 7 36 969];

end
%%%and adjust ewdev here
if beta>=0.93
ewdev(1:grid(3),:,:,3)=ewdev(index_s,:,:,3);
ewdev(grid(3)+1:end,:,:,3)=0;
wdev(1:grid(3),:,:,3)=wdev(index_s,:,:,3);
wdev(grid(3)+1:end,:,:,3)=0;
gamma3=gamma3_s;
end
gammasub(1:length(gamma3),1:3,3)=gamma3;
gammasub(1:length(gamma2),1:2,2)=gamma2;
%%%%and adjust ewdev for the penalty
for agentsub=2:3
for i=1:grid(agentsub)
for j=1:states(agentsub)
    for g=1:agentsub
    if ewdev(i,j,g,agentsub)>0
ewdev(i,j,g,agentsub)=ewdev(i,j,g,agentsub)*(1-penalty);
wdev(i,j,g,agentsub)=wdev(i,j,g,agentsub)*(1-penalty);
    else
ewdev(i,j,g,agentsub)=ewdev(i,j,g,agentsub)/(1-penalty);
wdev(i,j,g,agentsub)=wdev(i,j,g,agentsub)/(1-penalty);
    end
    end
end
end
end



gap=1;
gapcritsetexact=0.1
gapend=0.001;
options=optimset('Display','iter');
scalefactor=1;

count=0;
eps=0.0000001;
t=0;

wsbnew=[];
wsbnewsol=[];
phisol=[];
csol=[];
xkeep=[];
outputkeep=[];


tfinish=20
if beta==0.88
    tfinish=15
    
end
for t=1:tfinish

gap

agentperm=perms(1:agent);
count=0
countgrid=ones(n,state);
howmanybind=ones(n,state)*99;
countperm=ones(state,1);

X0=ones(4,2*agent);

%%%%%here we calculate the results when three agents have a binding
%%%%%constraint%%%%

h=agent-1;

countperm=ones(grid(agent-1),state);
countpermstate=ones(state,1);
for j=1:state
%%%%allocate binding and non-binding%%%%%
output(6)=0;
k=1;
bind=[1 2 3];
nobind=4;
%%%%allocate initial guesses%%%%%%%%%%%%%
for d=1:grid(agent-1)
d
if countperm(d,j)>0
output(6)=0;
%%%find a good starting guess

rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<=wdev(d,statedev(j,k,h),1,h) ...
    & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<=wdev(d,statedev(j,k,h),2,h) ...
    & utility(rho,cfirstbest(:,j,3))+beta*ewsb(:,j,3)<=wdev(d,statedev(j,k,h),3,h));
if size(rr)>0
x0=[cfirstbest(rr(end),j,1) cfirstbest(rr(end),j,2) cfirstbest(rr(end),j,3)];
x4=Y(j)-sum(x0); 
index=[1 j d k h];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);
if output(6)<=0
rrindex=1;
while rrindex<length(rr)   & output(6)<=0 
x0=[cfirstbest(rr(rrindex),j,1) cfirstbest(rr(rrindex),j,2) cfirstbest(rr(rrindex),j,3)];
index=[1 j d k h];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);
rrindex=rrindex+1;
end
end
end

%everything below are just different starting guesses
if output(6)<=0
   rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=wdev(d,statedev(j,k,h),1,h) ...
    & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=wdev(d,statedev(j,k,h),2,h) ...
    & utility(rho,cfirstbest(:,j,3))+beta*ewsb(:,j,3)>=wdev(d,statedev(j,k,h),3,h));
end
if size(rr)>0
x0=[cfirstbest(rr(1),j,1) cfirstbest(rr(1),j,2) cfirstbest(rr(1),j,3)];


for r=1:4
if output(6)<=0
index=[1 j d k h];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);
end
end
if output(6)<=0
rrindex=1;
while rrindex<length(rr)   & output(6)<=0 
x0=[cfirstbest(rr(rrindex),j,1) cfirstbest(rr(rrindex),j,2) cfirstbest(rr(rrindex),j,3)];
index=[1 j d k h];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);
rrindex=rrindex+1;
end
end
end
if output(6)<=0 & d>1
x0=[xkeep(d-1,j,1,k,h) xkeep(d-1,j,2,k,h) xkeep(d-1,j,3,k,h)   ];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);
end
if output(6)<=0
error('Optimal deviation did not converge');
end

outputkeepsol(d,j,1:agent,k,h)=output(1:4);
xkeepsol(d,j,1:agent,k,h)=[x(1:3) Y(j)-x(1)-x(2)-x(3)];
countperm(d,j)=0;
outputkeep(d,j,1:agent,k,h)=output(1:4);
xkeep(d,j,1:agent,k,h)=xkeepsol(d,j,1:agent,k,h);

%%%%allocate values within 1,2,3 binding

ipermod(1)=find(gamma3(:,1)==gamma3(d,1) & gamma3(:,2)==gamma3(d,3) & gamma3(:,3)==gamma3(d,2));
ipermod(2)=find(gamma3(:,1)==gamma3(d,2) & gamma3(:,2)==gamma3(d,1) & gamma3(:,3)==gamma3(d,3));
ipermod(3)=find(gamma3(:,1)==gamma3(d,2) & gamma3(:,2)==gamma3(d,3) & gamma3(:,3)==gamma3(d,1));
ipermod(4)=find(gamma3(:,1)==gamma3(d,3) & gamma3(:,2)==gamma3(d,1) & gamma3(:,3)==gamma3(d,2));
ipermod(5)=find(gamma3(:,1)==gamma3(d,3) & gamma3(:,2)==gamma3(d,2) & gamma3(:,3)==gamma3(d,1));

jpermod(1)=find(S(:,1)==S(j,1) & S(:,2)==S(j,3) & S(:,3)==S(j,2) & S(:,4)==S(j,4));
jpermod(2)=find(S(:,1)==S(j,2) & S(:,2)==S(j,1) & S(:,3)==S(j,3) & S(:,4)==S(j,4));
jpermod(3)=find(S(:,1)==S(j,2) & S(:,2)==S(j,3) & S(:,3)==S(j,1) & S(:,4)==S(j,4));
jpermod(4)=find(S(:,1)==S(j,3) & S(:,2)==S(j,1) & S(:,3)==S(j,2) & S(:,4)==S(j,4));
jpermod(5)=find(S(:,1)==S(j,3) & S(:,2)==S(j,2) & S(:,3)==S(j,1) & S(:,4)==S(j,4));

agentpermod(1,:)=[1 3 2 4];
agentpermod(2,:)=[2 1 3 4];
agentpermod(3,:)=[2 3 1 4];
agentpermod(4,:)=[3 1 2 4];
agentpermod(5,:)=[3 2 1 4];

for r=1:length(ipermod)
    for g=1:agent
    xkeep(ipermod(r),jpermod(r),g,k,h)=xkeepsol(d,j,agentpermod(r,g),k,h);
    outputkeep(ipermod(r),jpermod(r),g,k,h)=outputkeepsol(d,j,agentpermod(r,g),k,h);
    end
countperm(ipermod(r),jpermod(r))=0;
end


end  %end of permutations grid

end  %end of gamma grid

%%%binding for
%%%1 2 4,
%%%1 3 4,
%%%2 3 4
jperm_nck(1)=find(S(:,1)==S(j,1) & S(:,2)==S(j,2) & S(:,3)==S(j,4) & S(:,4)==S(j,3));
agentperm_nck(1,:)=[1 2 4 3];

jperm_nck(2)=find(S(:,1)==S(j,1) & S(:,2)==S(j,4) & S(:,3)==S(j,2) & S(:,4)==S(j,3));
agentperm_nck(2,:)=[1 4 2 3];

jperm_nck(3)=find(S(:,1)==S(j,4) & S(:,2)==S(j,1) & S(:,3)==S(j,2) & S(:,4)==S(j,3));
agentperm_nck(3,:)=[4 1 2 3];


for r=1:length(jperm_nck)
for g=1:agent
     outputkeep(:,jperm_nck(r),g,r+1,h)=outputkeep(:,j,agentperm_nck(r,g),1,h);
     xkeep(:,jperm_nck(r) ,g,r+1,h)=xkeep(:,j,agentperm_nck(r,g),1,h);
end
end


end   %%%end of state

clear outputkeepsol xkeepsol


%for j=1:state
tic


for j=1:state
stateperm=perms(S(j,:));
count=count+1

%%%%%now you need to check that the solutions are admissible, i.e. that no
%%%%%other coalition wants to deviate from it. 
%find the admissible deviations
%%%%consider deviations of 1,2 and 3, we want to check if combining with 4
%%%%gives a better deviation. Now you need check 1 2 4, 1 3 4, 2 3 4
%%%%i.e. keep 1,2 same, swap in 4 for 3; keep 1,3 same, swap in 4 for 2
%%%%etc.

admissible=[];
for d=1:grid(3)
    A=ones(grid(3),1)*[outputkeep(d,j,1,1,3) outputkeep(d,j,2,1,3)  outputkeep(d,j,4,1,3)];
    B=[outputkeep(1:grid(3),j,1,2,3) outputkeep(1:grid(3),j,2,2,3) outputkeep(1:grid(3),j,4,2,3)];
    C=sum(A+eps<=B,2);
    rr124=find(C<3);
    A=ones(grid(3),1)*[outputkeep(d,j,1,1,3) outputkeep(d,j,3,1,3)  outputkeep(d,j,4,1,3)];
    B=[outputkeep(1:grid(3),j,1,3,3) outputkeep(1:grid(3),j,3,3,3) outputkeep(1:grid(3),j,4,3,3)];
    C=sum(A+eps<=B,2);
    rr134=find(C<3);
    A=ones(grid(3),1)*[outputkeep(d,j,2,1,3) outputkeep(d,j,3,1,3)  outputkeep(d,j,4,1,3)];
    B=[outputkeep(1:grid(3),j,2,4,3) outputkeep(1:grid(3),j,3,4,3) outputkeep(1:grid(3),j,4,4,3)];
    C=sum(A+eps<=B,2);
    rr234=find(C<3);
   if size(rr124,1)==grid(3) & size(rr134,1)==grid(3) & size(rr234,1)==grid(3)
       admissible=[admissible; d];
   end
end
%%%%Now you need to do the same thing for admissible124 admissible134

admissible124=[];
for d=1:grid(3)
    A=ones(grid(3),1)*[outputkeep(d,j,1,2,3) outputkeep(d,j,2,2,3)  outputkeep(d,j,3,2,3)];
    B=[outputkeep(1:grid(3),j,1,1,3) outputkeep(1:grid(3),j,2,1,3) outputkeep(1:grid(3),j,3,1,3)];  %%%against 123
    C=sum(A+eps<=B,2);
    rr123=find(C<3);
    A=ones(grid(3),1)*[outputkeep(d,j,1,2,3) outputkeep(d,j,3,2,3)  outputkeep(d,j,4,2,3)];
    B=[outputkeep(1:grid(3),j,1,3,3) outputkeep(1:grid(3),j,3,3,3) outputkeep(1:grid(3),j,4,3,3)]; %%%against 134
    C=sum(A+eps<=B,2);
    rr134=find(C<3);
    A=ones(grid(3),1)*[outputkeep(d,j,2,2,3) outputkeep(d,j,3,2,3)  outputkeep(d,j,4,2,3)];
    B=[outputkeep(1:grid(3),j,2,4,3) outputkeep(1:grid(3),j,3,4,3) outputkeep(1:grid(3),j,4,4,3)];
    C=sum(A+eps<=B,2);
    rr234=find(C<3);
   if size(rr123,1)==grid(3) & size(rr134,1)==grid(3) & size(rr234,1)==grid(3)
       admissible124=[admissible124; d];
   end
end

%%%134
admissible134=[];
for d=1:grid(3)
    A=ones(grid(3),1)*[outputkeep(d,j,1,3,3) outputkeep(d,j,2,3,3)  outputkeep(d,j,3,3,3)];
    B=[outputkeep(1:grid(3),j,1,1,3) outputkeep(1:grid(3),j,2,1,3) outputkeep(1:grid(3),j,3,1,3)];  %%%against 123
    C=sum(A+eps<=B,2);
    rr123=find(C<3);
    A=ones(grid(3),1)*[outputkeep(d,j,1,3,3) outputkeep(d,j,2,3,3)  outputkeep(d,j,4,3,3)];
    B=[outputkeep(1:grid(3),j,1,2,3) outputkeep(1:grid(3),j,2,2,3) outputkeep(1:grid(3),j,4,2,3)]; %%%against 124
    C=sum(A+eps<=B,2);
    rr124=find(C<3);
    A=ones(grid(3),1)*[outputkeep(d,j,2,3,3) outputkeep(d,j,3,3,3)  outputkeep(d,j,4,3,3)];
    B=[outputkeep(1:grid(3),j,2,4,3) outputkeep(1:grid(3),j,3,4,3) outputkeep(1:grid(3),j,4,4,3)];
    C=sum(A+eps<=B,2);
    rr234=find(C<3);
   if size(rr123,1)==grid(3) & size(rr124,1)==grid(3) & size(rr234,1)==grid(3)
       admissible134=[admissible134; d];
   end
end





for i=1:n
 
index=[i,j]
if countgrid(i,j)>0

%%%%find which coalitions have binding constraints
for h=2:agent-1
TEMP=nchoosek(1:agent,h);
for k=1:size(TEMP,1)
dev(1:grid(h),k,h)=sum(ones(grid(h),1)*[squeeze(utility(rho,cfirstbest(i,j,TEMP(k,:)))+beta*ewsb(i,j,TEMP(k,:)))'] ....
    - [squeeze(wdev(1:grid(h),statedev(j,k,h),1:h,h))]<0,2);
if max(dev(1:grid(h),k,h))==h
BIND(k,1:h,h)=TEMP(k,1:h);
else
BIND(k,1:h,h)=0;
end
end
end
%%%%and find if there are any individual profitable deviations 
for g=1:agent
if (utility(rho,cfirstbest(i,j,g))+beta*ewsb(i,j,g))<(utility(rho,S(j,g))+beta*P(j,:)*waut(:,g))
BIND(g,1,1)=g;
else
BIND(g,1,1)=0;
end
end
    


%%%%find joint deviations
dev123=sum(ones(grid(3),1)*[utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1) utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2) utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3)]...
    -[outputkeep(1:grid(3),j,1,1,3) outputkeep(:,j,2,1,3) outputkeep(:,j,3,1,3)]<0,2);
dev124=sum(ones(grid(3),1)*[utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1) utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2) utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4)]...
    -[outputkeep(:,j,1,2,3) outputkeep(:,j,2,2,3) outputkeep(:,j,4,2,3)]<0,2);
dev134=sum(ones(grid(3),1)*[utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1) utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3) utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4)]...
    -[outputkeep(:,j,1,3,3) outputkeep(:,j,3,3,3) outputkeep(:,j,4,3,3)]<0,2);
dev234=sum(ones(grid(3),1)*[utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2) utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3) utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4)]...
    -[outputkeep(:,j,2,4,3) outputkeep(:,j,3,4,3) outputkeep(:,j,4,4,3)]<0,2);


coalitionbind=find(dev123==3);
%%%%and finding a binding constraint for two people

dev12=sum(ones(grid(2),1)*[utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1) utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)] ...
    - [ wdev(1:grid(2),statedev(j,1,2),1,2) ...
        wdev(1:grid(2),statedev(j,1,2),2,2)]<0,2);

dev13=sum(ones(grid(2),1)*[utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1) utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3)] ...
    - [ wdev(1:grid(2),statedev(j,2,2),1,2) ...
        wdev(1:grid(2),statedev(j,2,2),2,2)]<0,2);                        
                        
dev14=sum(ones(grid(2),1)*[utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1) utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4)] ...
    - [ wdev(1:grid(2),statedev(j,3,2),1,2) ...
        wdev(1:grid(2),statedev(j,3,2),2,2)]<0,2);                        
   
dev23=sum(ones(grid(2),1)*[utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2) utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3)] ...
    - [ wdev(1:grid(2),statedev(j,4,2),1,2) ...
        wdev(1:grid(2),statedev(j,4,2),2,2)]<0,2);                        
dev24=sum(ones(grid(2),1)*[utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2) utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4)] ...
    - [wdev(1:grid(2),statedev(j,5,2),1,2) ...
        wdev(1:grid(2),statedev(j,5,2),2,2)]<0,2);                        

dev34=sum(ones(grid(2),1)*[utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3) utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4)] ...
    - [ wdev(1:grid(2),statedev(j,6,2),1,2) ...
        wdev(1:grid(2),statedev(j,6,2),2,2)]<0,2);                        

%%%%this makes some adjustments that are needed because of the
%%%%interpolation
if S(j,2)==S(j,3) & gammagrid(i,2)+eps>=gammagrid(i,3) & gammagrid(i,2)-eps<=gammagrid(i,3)
    dev13=dev12;
end
if S(j,2)==S(j,4) & gammagrid(i,2)+eps>=gammagrid(i,4) & gammagrid(i,2)-eps<=gammagrid(i,4)
    dev14=dev12;
end
if S(j,1)==S(j,2) & gammagrid(i,1)==gammagrid(i,2)  & S(j,2)==S(j,3) & gammagrid(i,2)==gammagrid(i,3)
    dev23=dev12;
end
if S(j,1)==S(j,2) & gammagrid(i,1)==gammagrid(i,2)  & S(j,2)==S(j,4) & gammagrid(i,2)==gammagrid(i,4)
    dev24=dev12;
end
if S(j,1)==S(j,3) & gammagrid(i,1)==gammagrid(i,3)  & S(j,2)==S(j,4) & gammagrid(i,2)==gammagrid(i,4)
    dev34=dev12;
end


coalitionbind12=find(dev12==2);

%%this is specifically for t=8 & beta=0.96, j=13, i =177
if S(j,3)==S(j,4) & gammagrid(i,3)+eps>=gammagrid(i,4) & gammagrid(i,3)-eps<=gammagrid(i,4)  & max(dev12)==2
dev124(:)=2
end


   
gridperm=perms(gammagrid(i,1:4));
for k=1:length(gridperm)
iperm(k)=find(gammagrid(:,1)<=gridperm(k,1)+eps & gammagrid(:,1)>=gridperm(k,1)-eps& gammagrid(:,2)<=gridperm(k,2)+eps & gammagrid(:,2)>=gridperm(k,2)-eps & gammagrid(:,3)<=gridperm(k,3)+eps & gammagrid(:,3)>=gridperm(k,3)-eps & gammagrid(:,4)<=gridperm(k,4)+eps & gammagrid(:,4)>=gridperm(k,4)-eps);
jperm(k)=find(S(:,1)==stateperm(k,1) & S(:,2)==stateperm(k,2) & S(:,3)==stateperm(k,3) & S(:,4)==stateperm(k,4));
end

% Constraint binding for
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%   Person 1 only  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if scalefactor*(utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1))<scalefactor*(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2))>=scalefactor*(utility(rho,S(j,2))+beta*P(j,:)*waut(:,2)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))>=scalefactor*(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
    && scalefactor*(utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4))>=scalefactor*(utility(rho,S(j,4))+beta*P(j,:)*waut(:,4))...
    && max(dev12)<2  && max(dev13)<2 && max(dev14)<2 && max(dev23)<2 && max(dev24)<2 && max(dev34)<2 ...
    && max(dev123)<3  && max(dev124)<3 && max(dev134)<3 & max(dev234)<3;
bind=1;
nobind=[2 3 4];
bindsfor1(i,j)=1;
output(6)=0;

%%%would this even be possible? ex-ante checking
XBIND=[min(xkeep(1:grid(3),j,1,1,3)) min(xkeep(1:grid(3),j,2,1,3)) min(xkeep(1:grid(3),j,3,1,3)) min(xkeep(1:grid(3),j,4,2,3)) ] ;  
%%%%the first entry here is unimportant, second, 
%%%%third and fourth entry will be used to check if xaux (i.e. the likely
%%%%solution is less than what 2, 3 and 4 would get as a minimum in a
%%%%coalitional deviation
%%%The reason for doing these ex-ante checks is that otherwise you need to
%%%solve a problem that possibly goes off the grid but that would never be
%%%implemented anyway. 
xaux(1) = XBIND(1);
xaux(2) = (Y(j)-XBIND(1))/(1+gammagrid(i,7)/gammagrid(i,6)+gammagrid(i,8)/gammagrid(i,6));
xaux(3) = gammagrid(i,7)/gammagrid(i,6)*xaux(2);
xaux(4) = gammagrid(i,8)/gammagrid(i,6)*xaux(2);


bindpost=bind;

nobindpost=[];
for g=1:length(nobind)
     if xaux(nobind(g))<XBIND(nobind(g))
     bindpost=[bindpost nobind(g)];
     %xbindpost=[xbindpost XBIND(g)];
     else
     nobindpost=[nobindpost nobind(g)];
     end
end


%%%%%now go through the possibilities    
if length(bind)==length(bindpost)

index=[i j 1 1 1];
x0=[xaux(1:3)];
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
if output(6)<=0 & i>1
x0=[c(i-1,j,1) c(i-1,j,2) c(i-1,j,3)]
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
end
if output(6)<=0 
x0=[cfirstbest(i,j,1) cfirstbest(i,j,2) cfirstbest(i,j,3)]
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
end
if output(6)<=0
   error('agentsolverho_4agents did not converge') 
end
%%%check whether sharing the increase of 1 among 3 people there is a
%%%coalition of 3 or a coalition of 2 that would like to deviate

dev123_1=sum(ones(grid(3),1)*[output(1) output(2) output(3)] ...
       -squeeze(wdev(1:grid(3),statedev(j,1,3),1:3,3))<0,2);
dev124_1=sum(ones(grid(3),1)*[output(1) output(2) output(4)] ...
       -squeeze(wdev(1:grid(3),statedev(j,2,3),1:3,3))<0,2);
dev134_1=sum(ones(grid(3),1)*[output(1) output(3) output(4)] ...
       -squeeze(wdev(1:grid(3),statedev(j,3,3),1:3,3))<0,2);

%%%%and finding a binding constraint for two people
dev12_1=sum(ones(grid(2),1)*[output(1) output(2)]-squeeze(wdev(1:grid(2),statedev(j,1,2),1:2,2))<0,2);
dev13_1=sum(ones(grid(2),1)*[output(1) output(3)]-squeeze(wdev(1:grid(2),statedev(j,2,2),1:2,2))<0,2);
dev14_1=sum(ones(grid(2),1)*[output(1) output(4)]-squeeze(wdev(1:grid(2),statedev(j,3,2),1:2,2))<0,2);


    if max(dev12_1)<2  && max(dev13_1)<2 && max(dev14_1)<2  ...
    && ( size(dev123_1,1)==0 || max(dev123_1)<3 )  && (size(dev124_1,1)==0 || max(dev124_1)<3) && (size(dev134_1,1)==0 || max(dev134_1)<3 ); %%% noone else has a binding constraint
    bind=[1];
    nobind=[2 3 4];
    x=x(1:3);
    x(4)=Y(j)-x(1)-x(2)-x(3);
    
   
    %%%%first check that sharing person 1's increase in consumption it is not
    %%%%the case that 2 of the other 3 people have a binding constraint                       
    elseif max(dev123_1)==3  %%%person 1,2, and 3 have a binding constraint
    bind=[1 2 3];
    nobind=[4];
        if length(admissible)==1
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3)));
        else   
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3))');
        end
    x=squeeze(xkeep(admissible(d),j,1:4,1,3))';
    output=squeeze(outputkeep(admissible(d),j,1:4,1,3))';
    output(6)=1;
    
    elseif max(dev124_1)==3  %%%person 1,2, and 4 have a binding constraint
    bind=[1 2 4];
    nobind=[3];
        if length(admissible124)==1
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,2,3)));
        else   
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,2,3))');
        end
    x=squeeze(xkeep(admissible124(d),j,1:4,2,3))';
    output=squeeze(outputkeep(admissible124(d),j,1:4,2,3))';
    output(6)=1;
    elseif max(dev134_1)==3  %%%person 1,3, and 4 have a binding constraint
    bind=[1 3 4];
    nobind=[2];
        if length(admissible134)==1
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible134,j,:,3,3)));
        else   
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible134,j,:,3,3))');
        end
    x=squeeze(xkeep(admissible134(d),j,1:4,3,3))';
    output=squeeze(outputkeep(admissible134(d),j,1:4,3,3))';
    output(6)=1;
    %%%%now checking for any two: 
    elseif (max(dev12_1)==2 | max(dev13_1)==2 | max(dev14_1)==2)  ...
            && ( size(dev123_1,1)==0 || max(dev123_1)<3 )  && (size(dev124_1,1)==0 || max(dev124_1)<3) && (size(dev134_1,1)==0 || max(dev134_1)<3 );
        bindsfor12expost(i,j)=1;
        if max(dev12_1)==2   
        bind=[1 2];
        nobind=[3 4];
        k=1;
        elseif  max(dev13_1)==2 
        bind=[1 3];
        nobind=[2 4];
        k=2;
        elseif  max(dev14_1)==2
        bind=[1 4];
        nobind=[2 3];
        k=3;
        end
    %find all possible joint deviations of two people
    for d=1:grid(2)
        output(6)=0;

        index=[i j d k 2];
        x0=xaux;
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);
        %below are all different starting values
        if output(6)<=0 & d>1
        x0=[xkeepsol2(d-1,j,1) xkeepsol2(d-1,j,2) xkeepsol2(d-1,j,3)] ;
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        end
        if j==16 & output(6)<=0 
           x0=[XKEEPSOL2(d,1,1,i,k)  XKEEPSOL2(d,1,2,i,k) XKEEPSOL2(d,1,3,i,k)]*Y(16)/Y(1);
           [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub); 
        else
           [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);

        end
        if output(6)<=0 
        x0=[x];
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        end
        if output(6)<=0 && i>1 && size(XKEEPSOL2,2)>=j && size(XKEEPSOL2,5)>=k && XKEEPSOL2(d,j,1,i-1,k)~=0 
        x0=[XKEEPSOL2(d,j,1,i-1,k) XKEEPSOL2(d,j,2,i-1,k) XKEEPSOL2(d,j,3,i-1,k) ];
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        end
        
        if output(6)<=0 
        nobindaux=[];
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindaux,gammasub);
        
        x0(1:2)=x(1:2);
        x0(3) = (Y(j)-x(1)-x(2))/(1+gammagrid(i,8)/gammagrid(i,7));
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        end
        if output(6)<=0
        if j==8
        x0=[1.0154 1.4514 1.8976];
        end
        [x, output] = maxall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        x0=x;
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        end

        if output(6)<=0
            for kk=1:5
                tol=10^(-(kk+1))
               [x, output] = solveall_tol_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,tol,gammasub);
                x0=x;
            exitflag(d,j,kk,i)=output(6);
            exitflag(d,j,kk,i)
            end
            if exitflag(d,j,1,i)<=0
            output(6)=0
        else
            output(6)=1
            end

        end
        
        
        outputkeepsol2(d,j,1:agent)=output(1:4);
        xkeepsol2(d,j,1:agent)=[x(1:3) Y(j)-x(1)-x(2)-x(3)];
        XKEEPSOL2(d,j,1:agent,i,k)=[x(1:3) Y(j)-x(1)-x(2)-x(3)];
      
        if output(6)<=0
        error('did not converge for 1 and 2 binding')
        end
    end
        %find the admissible ones
    admissible2=[]
            for d=1:grid(2)
            dev123_2=sum(ones(grid(3),1)*[outputkeepsol2(d,j,1) outputkeepsol2(d,j,2) outputkeepsol2(d,j,3)]...
                -squeeze(wdev(1:grid(3),statedev(j,1,3),1:3,3))<0,2);
            dev124_2=sum(ones(grid(3),1)*[outputkeepsol2(d,j,1) outputkeepsol2(d,j,2) outputkeepsol2(d,j,4)]...
                -squeeze(wdev(1:grid(3),statedev(j,2,3),1:3,3))<0,2);


          
            r123=find(dev123_2==3);
            r124=find(dev124_2==3);
            if size(r123)>0 | size(r124)>0
            else
            admissible2=[admissible2;d]; 
            end
            end

        %%%%Step 2: find the optimal ones (if none are admissible, need to go to
        %%%%solution for 3 binding)

            if size(admissible2)>0 
                if size(admissible2,1)>1 
                [CC,d]=max([(gammagrid(i,1)) gammagrid(i,2:4)]*squeeze(outputkeepsol2(admissible2,j,:))');
                else
                [CC,d]=max([(gammagrid(i,1)) gammagrid(i,2:4)]*squeeze(outputkeepsol2(admissible2,j,:)));
                end
            x=xkeepsol2(admissible2(d),j,:);
            output=outputkeepsol2(admissible2(d),j,:);
            output(6)=1;
            elseif size(r123)>0
            bind=[1 2 3];
            nobind=4;
            if length(admissible)==1
            [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3)));
            else   
            [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3))');
            end
            x=xkeep(admissible(d),j,:,1,3);
            output=outputkeep(admissible(d),j,:,1,3);
            output(6)=1;
            elseif size(r124)>0
            bind=[1 2 4];
            nobind=3;
            if length(admissible124)==1
            [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,2,3)));
            else   
            [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,2,3))');
            end
            x=xkeep(admissible124(d),j,:,2,3);
            output=outputkeep(admissible124(d),j,:,2,3);
            output(6)=1;    
            end
    elseif max(dev12_1)<2  && max(dev13_1)==2 && max(dev14_1)<2  ...
    && max(dev123_1)<3  && max(dev124_1)<3 && max(dev134_1)<3;
        error('permutation is not working')
    elseif max(dev12_1)<2  && max(dev13_1)<2 && max(dev14_1)==2  ...
    && max(dev123_1)<3  && max(dev124_1)<3 && max(dev134_1)<3; %%
        error('permutation is not working')
    end
elseif length(bind)<length(bindpost)
    %%%%now checking for any two who are constrained 
       
    if length(bindpost)==2
        bindsfor1bindpost2(i,j)=1;
        
        if bindpost==[1 2]
            k=1;
        elseif bindpost==[1 3]
            k=2;
        elseif bindpost==[1 4]
            k=3;    
        end
        for d=1:grid(2)
        output(6)=0;

        index=[i j d k 2];
        rr=find(utility(rho,cfirstbest(:,j,bindpost(1)))+beta*ewsb(:,j,bindpost(1))<=wdev(d,statedev(j,k,2),1,2) ...
        & utility(rho,cfirstbest(:,j,bindpost(2)))+beta*ewsb(:,j,bindpost(2))<=wdev(d,statedev(j,k,2),1,2) );
        if k==1
        x0(1)=[cfirstbest(rr(end),j,bindpost(1))];
        x0(2)=[cfirstbest(rr(end),j,bindpost(2))];
        x0(3) = (Y(j)-x0(bindpost(1))-x0(bindpost(2)))/(1+gammagrid(i,nobindpost(2)+4)/gammagrid(i,nobindpost(1)+4));
        elseif k==2
        x0(1)=[cfirstbest(rr(end),j,bindpost(1))];
        x0(3)=[cfirstbest(rr(end),j,bindpost(2))];
        x0(2) = (Y(j)-x0(bindpost(1))-x0(bindpost(2)))/(1+gammagrid(i,nobindpost(2)+4)/gammagrid(i,nobindpost(1)+4));
        elseif k==3
        x0(1)=[cfirstbest(rr(end),j,bindpost(1))];
        x0(2) = (Y(j)-x0(bindpost(1))-cfirstbest(rr(end),j,bindpost(2)))/(1+gammagrid(i,nobindpost(2)+4)/gammagrid(i,nobindpost(1)+4));
        x0(3)=[Y(j)-x0(1)-x0(2)-cfirstbest(rr(end),j,bindpost(2))];
        end
        
            
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        %below are all different starting values
        if output(6)<=0 & d>1
        x0=[xkeepsol2(d-1,j,1) xkeepsol2(d-1,j,2) xkeepsol2(d-1,j,3)] ;
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        end
        if j==16
           x0=[XKEEPSOL2(d,1,1,i,k)  XKEEPSOL2(d,1,2,i,k) XKEEPSOL2(d,1,3,i,k)]*Y(16)/Y(1);
           [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub); 
        else
           [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);

        end
        if output(6)<=0 
        x0=[x];
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        end
        if output(6)<=0 && i>1 && size(XKEEPSOL2,2)>=j && size(XKEEPSOL2,5)>=k && XKEEPSOL2(d,j,1,i-1,k)~=0 
        x0=[XKEEPSOL2(d,j,1,i-1,k) XKEEPSOL2(d,j,2,i-1,k) XKEEPSOL2(d,j,3,i-1,k) ];
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,gammasub);
        end
        
        if output(6)<=0 
        nobindaux=[];
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindaux, gammasub);
        
        x0(1:2)=x(1:2);
        x0(3) = (Y(j)-x(1)-x(2))/(1+gammagrid(i,8)/gammagrid(i,7));
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost, gammasub);
        end
        if output(6)<=0
        if j==8
        x0=[1.0154 1.4514 1.8976];
        end
        [x, output] = maxall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost, gammasub);
        x0=x;
        [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost, gammasub);
        end

        if output(6)<=0
            for kk=1:5
                tol=10^(-(kk+1))
               [x, output] = solveall_tol_shareimmediate(x0,index,ewsb,waut,bindpost,nobindpost,tol,gammasub);
                x0=x;
            exitflag(d,j,kk,i)=output(6);
            exitflag(d,j,kk,i)
            end
            if exitflag(d,j,1,i)<=0
            output(6)=0
            else
            output(6)=1
            end

        end
        
        outputkeepsol2(d,j,1:agent)=output(1:4);
        xkeepsol2(d,j,1:agent)=[x(1:3) Y(j)-x(1)-x(2)-x(3)];
        XKEEPSOL2(d,j,1:agent,i,k)=[x(1:3) Y(j)-x(1)-x(2)-x(3)];
        if output(6)<=0
        error('did not converge for 1 and 2 binding')
        end
        end
       
        
        
       %%%%Step 2: find the admissible ones
        admissible2=[]


            for d=1:grid(2)
            dev123_2=sum(ones(length(admissible),1)*[outputkeepsol2(d,j,1) outputkeepsol2(d,j,2) outputkeepsol2(d,j,3)]...
                -squeeze(wdev(1:grid(3),statedev(j,1,3),1:3,3))<0,2);
            dev124_2=sum(ones(length(admissible124),1)*[outputkeepsol2(d,j,1) outputkeepsol2(d,j,2) outputkeepsol2(d,j,4)]...
                -squeeze(wdev(1:grid(3),statedev(j,2,3),1:3,3))<0,2);
            
            r123=find(dev123_2==3)
            r124=find(dev124_2==3);
            if size(r124)>0 | size(r123)>0
            else
            admissible2=[admissible2;d]; 
            end
            end

        %%%%Step 2: find the optimal ones (if none are admissible, need to go to
        %%%%solution for 3 binding)

            if size(admissible2)>0 
                if size(admissible2,1)>1  
                [CC,d]=max([(gammagrid(i,1)) gammagrid(i,2:4)]*squeeze(outputkeepsol2(admissible2,j,:))');
                else
                [CC,d]=max([(gammagrid(i,1)) gammagrid(i,2:4)]*squeeze(outputkeepsol2(admissible2,j,:)));
                end
            x=xkeepsol2(admissible2(d),j,:);
            output=outputkeepsol2(admissible2(d),j,:);
            output(6)=1;
            elseif size(r123)>0
            bind=[1 2 3];
            nobind=4;
            if length(admissible)==1
            [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3)));
            else   
            [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3))');
            end
            x=squeeze(xkeep(admissible(d),j,:,1,3))';
            output=squeeze(outputkeep(admissible(d),j,:,1,3))';
            output(6)=1;
            elseif size(r124)>0
            bind=[1 2 4];
            nobind=3;
            if length(admissible124)==1
            [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,2,3)));
            else   
            [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,2,3))');
            end
            x=squeeze(xkeep(admissible124(d),j,:,2,3))';
            output=squeeze(outputkeep(admissible124(d),j,:,2,3))';
            output(6)=1;    
            end
        
    elseif length(bindpost)==3
        %%%%
        k=find(bindpost(1)==TEMP(:,1) & bindpost(2)==TEMP(:,2) & bindpost(3)==TEMP(:,3));
        if k<=3
            if k==1
        if length(admissible)==1
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,k,3)));
        else   
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,k,3))');
        end
        x=squeeze(xkeep(admissible(d),j,:,k,3))';
        output=squeeze(outputkeep(admissible(d),j,:,k,3))';
        output(6)=1; 
            elseif k==2
                if length(admissible124)==1
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,k,3)));
        else   
        [CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,k,3))');
        end
        x=squeeze(xkeep(admissible124(d),j,:,k,3))';
        output=squeeze(outputkeep(admissible124(d),j,:,k,3))';
        output(6)=1; 
            end
                

        elseif k>3
        error('Permutation does not work')
        end  
        
        
    end
    
end



csol(i,j,1:agent-1)=x(1:agent-1);
csol(i,j,agent)=Y(j)-sum(x(1:agent-1));
wsbnewsol(i,j,1:agent)=output(1:4);

for g=1:agent
gammar(i,j,g)=csol(i,j,g)/csol(i,j,1);
end
for g=1:agent
phisol(i,j,g)=((gammagrid(i,nobind(end))/gammagrid(i,g))/(gammar(i,j,nobind(end))/gammar(i,j,g)))^rho-1;
end

%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
howmanybind(iperm(k),jperm(k))=1;
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Constraint binding for Person 1 and Person 2 %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif (size(coalitionbind12,1)>0 && size(coalitionbind12,2)>0 ...  
    && scalefactor*(utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))>=scalefactor*(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
    && scalefactor*(utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4))>=scalefactor*(utility(rho,S(j,4))+beta*P(j,:)*waut(:,4))...
    && max(dev123)<3  && max(dev134)<3 && max(dev124)<3 && max(dev234)<3) ;   
  bind=[1 2];
  nobind=[3 4];

  bindsfor12(i,j)=1;

%%%%Ex-ante checking
XBIND=[min(xkeep(1:grid(3),j,1,1,3)) min(xkeep(1:grid(3),j,2,1,3)) min(xkeep(1:grid(3),j,3,1,3)) min(xkeep(1:grid(3),j,4,2,3)) ] ;  %%%%the first two entries here are unimportant, 
%%%%third and fourth entry will be used to check if xaux (i.e. the likely
%%%%solution is less than what 3 and 4 would get as a minimum in a coalitional deviation
%%%% The reason for doing these ex-ante checks is that otherwise Matlab
%%%% needs to solve a complicated problem that would never be implemented. 
xaux(1,1:2)=XBIND(1:2);
xaux(3) = (Y(j)-xaux(1)-xaux(2))/(1+gammagrid(i,8)/gammagrid(i,7));
xaux(4) = gammagrid(i,8)/gammagrid(i,7)*xaux(3);

bindpost=bind;

nobindpost=[];
for g=1:length(nobind)
     if xaux(nobind(g))<XBIND(nobind(g))
     bindpost=[bindpost nobind(g)];
    
     else
     nobindpost=[nobindpost nobind(g)];
     end
end


if length(bind)==length(bindpost) | size(admissible)==0
    bindequalbindpost12(i,j)=1;
    
    dgrid=[1:7];
%find all the possible divisions of surplus for 1 and 2 joint deviation
for d=dgrid
output(6)=0;
%%%find a good starting guess
rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<=wdev(d,statedev(j,1,2),1,2) ...
    & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<=wdev(d,statedev(j,1,2),2,2) );
x0=[cfirstbest(rr(end),j,1) cfirstbest(rr(end),j,2) cfirstbest(rr(end),j,3)];
x0(3) = (Y(j)-x0(1)-x0(2))/(1+gammagrid(i,8)/gammagrid(i,7));
index=[i j d 1 2];
 [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
%%this is all for finding different starting values
if output(6)<=0
rrindex=1;
while rrindex<length(rr) & output(6)<=0    
x0=[cfirstbest(rr(rrindex),j,1) cfirstbest(rr(rrindex),j,2) cfirstbest(rr(rrindex),j,3)];
x0(3) = (Y(j)-x0(1)-x0(2))/(1+gammagrid(i,8)/gammagrid(i,7));
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
rrindex=rrindex+1;
end
end
if j==16 & output(6)<=0 
   x0=[XKEEPSOL2(d,1,1,i)  XKEEPSOL2(d,1,2,i) XKEEPSOL2(d,1,3,i)]*Y(16)/Y(1);
   [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub); 
else
   [x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
end

if output(6)<=0 
x0=[x];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
end
if output(6)<=0 && i>1 && size(XKEEPSOL2,2)>=j && XKEEPSOL2(d,j,1,i-1)~=0 
x0=[XKEEPSOL2(d,j,1,i-1) XKEEPSOL2(d,j,2,i-1) XKEEPSOL2(d,j,3,i-1) ];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);
end
if output(6)<=0 & d>1
x0=[xkeepsol2(d-1,j,1) xkeepsol2(d-1,j,2) xkeepsol2(d-1,j,3)];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
end
if output(6)<=0 
nobind=[];
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
nobind=[3 4];
x0(1:2)=x(1:2);
x0(3) = (Y(j)-x(1)-x(2))/(1+gammagrid(i,8)/gammagrid(i,7));
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
end
if output(6)<=0
if j==8
x0=[1.0154 1.4514 1.8976];
end
[x, output] = maxall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);
x0=x;
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
end
if output(6)<=0 & t>1
x0=[c(i,j,1) c(i,j,2) c(i,j,3)]    
[x, output] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
end

if output(6)<=0
    for kk=1:5
        tol=10^(-(kk+1))
       [x, output] = solveall_tol_shareimmediate(x0,index,ewsb,waut,bind,nobind,tol, gammasub);
        x0=x;
    exitflag(d,j,kk,i)=output(6);
    exitflag(d,j,kk,i)
    end
    if exitflag(d,j,1,i)<=0 & beta==0.97
    [x, output] = maxall_shareimmediate(x0,index,ewsb,waut,bind,nobind,gammasub);

    output(6)=1
    else
    output(6)=1
    rr=find(exitflag(d,j,:,i)>0)
    tol=10^(-(rr(end)+1));
    [x, output] = solveall_tol_shareimmediate(x0,index,ewsb,waut,bind,nobind,tol,gammasub);
    
    [x1, output1] = maxall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
    x0=x1;
    [x2, output2] = solveall_shareimmediate(x0,index,ewsb,waut,bind,nobind, gammasub);
    if output2(6)>0
        x=x2;
        output=output2;
    end
    
    end
    
end


outputkeepsol2(d,j,1:agent)=output(1:4);
xkeepsol2(d,j,1:agent)=[x(1:3) Y(j)-x(1)-x(2)-x(3)];
XKEEPSOL2(d,j,1:agent,i,1)=[x(1:3) Y(j)-x(1)-x(2)-x(3)];
if output(6)<=0
error('did not converge for 1 and 2 binding')
end
end
%%Check which of those are admissible 
admissible2=[]

for d=dgrid
dev123_2=sum(ones(length(admissible),1)*[outputkeepsol2(d,j,1) outputkeepsol2(d,j,2) outputkeepsol2(d,j,3)]...
    -squeeze(wdev(1:grid(3),statedev(j,1,3),1:3,3))<0,2);
dev124_2=sum(ones(length(admissible124),1)*[outputkeepsol2(d,j,1) outputkeepsol2(d,j,2) outputkeepsol2(d,j,4)]...
    -squeeze(wdev(1:grid(3),statedev(j,2,3),1:3,3))<0,2);

r123=find(dev123_2==3);
r124=find(dev124_2==3);
if size(r123)>0 | size(r124)>0
else
admissible2=[admissible2;d]; 
end
end

%%%%Step 2: find the optimal ones (if none are admissible, need to go to
%%%%solution for 3 binding, i.e. either 123 or 124)

if size(admissible2)>0 
    if size(admissible2,1)>1  
    [CC,d]=max([(gammagrid(i,1)) gammagrid(i,2:4)]*squeeze(outputkeepsol2(admissible2,j,:))');
    else
    [CC,d]=max([(gammagrid(i,1)) gammagrid(i,2:4)]*squeeze(outputkeepsol2(admissible2,j,:)));
    end
x=squeeze(xkeepsol2(admissible2(d),j,:))';
output=squeeze(outputkeepsol2(admissible2(d),j,:))';
output(6)=1;
elseif size(r123)>0
bind=[1 2 3];
nobind=4;
if length(admissible)==1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3)));
else   
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3))');
end
x=squeeze(xkeep(admissible(d),j,:,1,3))';
output=squeeze(outputkeep(admissible(d),j,:,1,3))';
output(6)=1;
elseif size(r124)>0
bind=[1 2 4];
nobind=3;
[CC,d]=max([(gammagrid(i,1)) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,2,3))');

x=squeeze(xkeep(admissible124(d),j,:,2,3))';
output=squeeze(outputkeep(admissible124(d),j,:,2,3))';
output(6)=1;    
end    


else
k=find(bindpost(1)==TEMP(:,1) & bindpost(2)==TEMP(:,2) & bindpost(3)==TEMP(:,3));
if k==1
if length(admissible)==1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,k,3)));
else   
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,k,3))');
end
x=squeeze(xkeep(admissible(d),j,:,k,3))';
output=squeeze(outputkeep(admissible(d),j,:,k,3))';
elseif k==2
if length(admissible124)==1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,k,3)));
else   
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible124,j,:,k,3))');
end
x=squeeze(xkeep(admissible124(d),j,:,k,3))';
output=squeeze(outputkeep(admissible124(d),j,:,k,3))';   
elseif k==3
if length(admissible)==1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible134,j,:,k,3)));
else   
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible134,j,:,k,3))');
end
x=squeeze(xkeep(admissible134(d),j,:,k,3))';
output=squeeze(outputkeep(admissible134(d),j,:,k,3))';    
end
output(6)=1;    

end

if output(6)<=0
error('Did not converge');
end

%%%allocate the solutions    
csol(i,j,1:agent-1)=x(1:agent-1);
csol(i,j,agent)=Y(j)-sum(x(1:agent-1));
wsbnewsol(i,j,1:agent)=output(1:4);

for g=1:agent
gammar(i,j,g)=csol(i,j,g)/csol(i,j,1);
end
for g=1:agent
phisol(i,j,g)=((gammagrid(i,nobind(end))/gammagrid(i,g))/(gammar(i,j,nobind(end))/gammar(i,j,g)))^rho-1;
end


%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
howmanybind(iperm(k),jperm(k))=2;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constraint binding for Person 1 and Person 2 and Person 3%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif size(coalitionbind,1)>0 & size(coalitionbind,2)>0  ;
bind=[1 2 3];
nobind=4;
if length(admissible)==1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3)));
else  
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:4)]*squeeze(outputkeep(admissible,j,:,1,3))');
end
history(i,j)=admissible(d);
csol(i,j,1:4)=xkeep(admissible(d),j,1:4,1,3);
wsbnewsol(i,j,1:agent)=outputkeep(admissible(d),j,1:agent,1,3);
countgrid(i,j)=0;

for g=1:agent
gammar(i,j,g)=csol(i,j,g)/csol(i,j,1);
end
for g=1:agent
phisol(i,j,g)=((gammagrid(i,nobind(end))/gammagrid(i,g))/(gammar(i,j,nobind(end))/gammar(i,j,g)))^rho-1;
end

%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
howmanybind(iperm(k),jperm(k))=3;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%Constraint binding for noone%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif scalefactor*(utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1))>=scalefactor*(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2))>=scalefactor*(utility(rho,S(j,2))+beta*P(j,:)*waut(:,2)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))>=scalefactor*(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
    && scalefactor*(utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4))>=scalefactor*(utility(rho,S(j,4))+beta*P(j,:)*waut(:,4)) ...
    && max(dev12)<2  && max(dev13)<2 && max(dev14)<2 && max(dev23)<2 && max(dev24)<2 && max(dev34)<2 ...
    && max(dev123)<3  && max(dev124)<3 && max(dev134)<3 && max(dev234)<3;
bind=[];
nobind=[1:4];
csol(i,j,:)=cfirstbest(i,j,:);
for g=1:agent
wsbnewsol(i,j,g)=utility(rho,cfirstbest(i,j,g))+beta*ewsb(i,j,g);
end
countgrid(i,j)=0;

for g=1:agent
gammar(i,j,g)=csol(i,j,g)/csol(i,j,1);
end
for g=1:agent
phisol(i,j,g)=((gammagrid(i,nobind(end))/gammagrid(i,g))/(gammar(i,j,nobind(end))/gammar(i,j,g)))^rho-1;
end

%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
howmanybind(iperm(k),jperm(k))=0;
end


end




end % end of countgrid if loop
end  % end of i loop


end  % end of j loop
toc
gap=scalefactor*norm((wsbnew(:,:,1)-wsb(:,:,1))./wsb(:,:,1));
gap
%gamma's
for j=1:state
for g=1:agent
gammar(:,j,g)=c(:,j,g)./c(:,j,1);
end
end
if sum(sum(countgrid))>0
    error('countgrid not zero')
end
wsb=wsbnew;
for j=1:state
        for g=1:agent
ewsb(:,j,g)=wsb(:,:,g)*P(j,:)';
        end
end
gap
GAP(t)=gap;
clear exitflag

end % of iteration loop



